---
title: "Provincial analysis"
author: "Francesco Bailo"
date: "2 July 2016"
output: pdf_document
---

```{r}
require(knitr)
setwd("~/Documents/replication_packages/eps_2016/")
opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE) 
```


```{r}
immi <- 
  read.csv('provincial_immigration_data.csv',
           stringsAsFactors = FALSE)
votes <- 
  read.csv('provincial_voting_data.csv',
           stringsAsFactors = FALSE)
votes$year <- gsub('\\d{2}_\\d{2}_', "", votes$data) 

# Eco stats
## Import province
codici_province <- 
  read.table("istat_province_code.csv", 
             sep=";", quote="\"", stringsAsFactors=FALSE, header = TRUE, encoding = 'latin1')
codici_province <-
  codici_province[,c('Codice.Regione','Codice.Provincia..1.',
                     'Denominazione.Città.metropolitana', 
                     'Denominazione.provincia')]
codici_province[[3]][codici_province[[3]] == "-"] <- NA
codici_province[[4]][codici_province[[4]] == "-"] <- NA
col_names <- c('cod_regione', 'cod_provincia', 'nome_provincia')
df1 <- subset(codici_province[!is.na(codici_province[[3]]),], 
              select = c('Codice.Regione','Codice.Provincia..1.',
                         'Denominazione.Città.metropolitana'))
names(df1) <- col_names
df1 <- unique(df1)

df2 <- subset(codici_province[!is.na(codici_province[[4]]),], 
              select = c('Codice.Regione','Codice.Provincia..1.',
                         'Denominazione.provincia'))
names(df2) <- col_names
df2 <- unique(df2)

codici_province <- rbind(df1, df2)
codici_province <- 
  codici_province[complete.cases(codici_province),]
codici_province$nome_provincia <- 
  gsub("/(.*)", "", codici_province$nome_provincia)

codici_province$nome_provincia <- toupper(codici_province$nome_provincia)

## Unemployment
unempl <-  
  read.csv("istat_province_unempl.csv", 
           encoding = 'latin1', comment.char="#", stringsAsFactors = FALSE)
unempl <- subset(unempl, 
                 Classe.di.età == '15 anni e più' & Sesso == 'totale')
unempl$Territorio[unempl$Territorio == "Valle d'Aosta / Vallée d'Aoste"] <- 
  "Valle d'Aosta"
unempl$Territorio[unempl$Territorio == "Provincia Autonoma Bolzano / Bozen"] <- 
  "Bolzano"
unempl$Territorio[unempl$Territorio == "Provincia Autonoma Trento"] <- 
  "Trento"
unempl$Territorio <- toupper(unempl$Territorio)

unempl$Territorio[!(unempl$Territorio %in% codici_province$nome_provincia)]

unempl <- merge(unempl, codici_province, 
                by.x = 'Territorio', 
                by.y = 'nome_provincia')
unempl <- unempl[,c('Tempo.e.frequenza', 'X0', 'cod_provincia')]
names(unempl) <- c('year', 'perc_unempl', 'cod_provincia')

# a2 <- data.frame(reg = as.character(immi$code_region),
#                  prov = as.character(immi$code_province),
#                  year = as.character(immi$year))
# a1 <- data.frame(reg = as.character(votes$cod_regione),
#                  prov = as.character(votes$cod_provincia),
#                  year = as.character(votes$year))
# 
# v1 <- apply(a1, 1, FUN = function(x) paste(x, collapse = '---'))
# v2 <- apply(a2, 1, FUN = function(x) paste(x, collapse = '---'))
# 
# v1[!(v1 %in% v2)]

# Province of Monza is created in 2009, 
# all data concerning Monza is summed to the province of Milan

immi$province[immi$code_province == 108] <- 'Milano'
immi$code_province[immi$code_province == 108] <- 15
require(dplyr)

# cat(paste(paste0(names(immi), " = sum(", names(immi),")"), collapse = ', \n'))

immi <- 
  immi %>%
  dplyr::group_by(province, year, code_region, code_province) %>%
  dplyr::summarize(foreign_residents = sum(foreign_residents),
                   male_foreign_residents = sum(male_foreign_residents),
                   female_foreign_residents = sum(female_foreign_residents),
                   muslim_foreign_residents = sum(muslim_foreign_residents),
                   Australia.and.New.Zealand = sum(Australia.and.New.Zealand), 
                   Caribbean = sum(Caribbean), 
                   Central.America = sum(Central.America), 
                   Central.Asia = sum(Central.Asia), 
                   Eastern.Africa = sum(Eastern.Africa), 
                   Eastern.Asia = sum(Eastern.Asia), 
                   Eastern.Europe = sum(Eastern.Europe), 
                   Melanesia = sum(Melanesia), 
                   Micronesia = sum(Micronesia), 
                   Middle.Africa = sum(Middle.Africa), 
                   Northern.Africa = sum(Northern.Africa), 
                   Northern.America = sum(Northern.America), 
                   Northern.Europe = sum(Northern.Europe), 
                   Polynesia = sum(Polynesia), 
                   South.America = sum(South.America), 
                   South.Eastern.Asia = sum(South.Eastern.Asia), 
                   Southern.Africa = sum(Southern.Africa), 
                   Southern.Asia = sum(Southern.Asia), 
                   Southern.Europe = sum(Southern.Europe), 
                   Western.Africa = sum(Western.Africa), 
                   Western.Asia = sum(Western.Asia), 
                   Western.Europe = sum(Western.Europe))

votes$Ente[votes$cod_provincia == 108] <- "MILANO"
votes$cod_provincia[votes$cod_provincia == 108] <- 15
votes <-
  votes %>%
  dplyr::group_by(Ente, cod_regione, cod_provincia, year, data, elezione) %>%
  dplyr::summarize(Numero.elettori = sum(Numero.elettori),
                   Numero.votanti = sum(Numero.votanti),
                   Schede.bianche = sum(Schede.bianche),
                   Schede.non.valide = sum(Schede.non.valide),
                   Voti.LN = sum(Voti.LN))

# Add age groups 
load("istat_pop_province_by_age.RData")
pop_by_age$over65Pop <- pop_by_age$over65_pop / pop_by_age$tot_pop
votes <- 
  merge(votes, pop_by_age[,c('cod_prov','year', 'over65Pop')], 
        all.x = TRUE, 
        by.x = c('cod_provincia','year'),
        by.y = c('cod_prov','year'))

immi_votes <-
  merge(immi, votes, 
        by.x = c('code_region','code_province','year'),
        by.y = c('cod_regione','cod_provincia','year'),
        all.x = TRUE, all.y = TRUE)

immi_votes$perc_pop_straniera <- 
  with(immi_votes, foreign_residents / Numero.elettori)
immi_votes$perc_pop_female <- 
  with(immi_votes, female_foreign_residents / Numero.elettori)
immi_votes$perc_pop_male <- 
  with(immi_votes, male_foreign_residents / Numero.elettori)
immi_votes$perc_pop_musul <- 
  with(immi_votes, muslim_foreign_residents / Numero.elettori)
immi_votes$perc_pop_non_musul <- 
  with(immi_votes, (foreign_residents - muslim_foreign_residents) / Numero.elettori)

# cat(paste0(names(immi)[7:28], collapse = '", "'))

region_vars = c("Australia.and.New.Zealand", "Caribbean", "Central.America", "Central.Asia",
                "Eastern.Africa", "Eastern.Asia", "Eastern.Europe", "Melanesia",
                "Micronesia", "Middle.Africa", "Northern.Africa", "Northern.America",
                "Northern.Europe", "Polynesia", "South.America", "South.Eastern.Asia",
                "Southern.Africa", "Southern.Asia", "Southern.Europe", "Western.Africa",
                "Western.Asia", "Western.Europe")

for(r in region_vars) {
  immi_votes[[paste0('perc_', r)]] <- immi_votes[[r]] / immi_votes$Numero.elettori
}

immi_votes$perc_lega <- 
  with(immi_votes, Voti.LN / (Numero.votanti - Schede.bianche - Schede.non.valide))
immi_votes$affluenza <- 
  with(immi_votes, Numero.votanti / Numero.elettori)

immi_votes <- unique(immi_votes)
immi_votes$code_province <- as.integer(immi_votes$code_province)

immi_votes <- 
  immi_votes %>%
  dplyr::group_by(code_region, year) %>%
  dplyr::mutate(perc_lega_reg = sum(Voti.LN) / sum((Numero.votanti - Schede.bianche - Schede.non.valide)))

immi_votes <- 
  merge(immi_votes, unempl, 
        by.x = c('year', 'code_province'),
        by.y = c('year', 'cod_provincia'), all.x = TRUE)


# Immigration trends
require(ggplot2)
immi_by_year <-
  immi %>%
  dplyr::group_by(year) %>%
  dplyr::summarize(foreign_residents  = sum(foreign_residents),
                   muslim_foreign_residents = sum(muslim_foreign_residents))

require(reshape2)
immi_by_year_by_world_region <-
  melt(immi, id.vars = c('year', 'province', 'code_region', 'code_province',
                         'foreign_residents', 'muslim_foreign_residents')) %>% 
  dplyr::group_by(year, variable) %>%
  dplyr::summarize(foreign_residents = sum(value))

ggplot() + 
  geom_line(data = immi_by_year, aes(x=year,y=foreign_residents), colour = 'blue') +
  geom_point(data = immi_by_year, aes(x=year,y=foreign_residents), colour = 'blue') +
  geom_line(data = immi_by_year, aes(x=year,y=muslim_foreign_residents), colour = 'red') +
  geom_point(data = immi_by_year, aes(x=year,y=muslim_foreign_residents), colour = 'red') +
  theme_bw() + labs(x='Foreign residents (blue); Muslim foreign residents (red)')

average_foreign_pop_world_regions <-
  immi_by_year_by_world_region %>%
  dplyr::group_by(variable) %>%
  dplyr::summarize(mean = mean(foreign_residents))
average_foreign_pop_world_regions$mean <- 
  with(average_foreign_pop_world_regions, mean / sum(mean))
average_foreign_pop_world_regions <-
  average_foreign_pop_world_regions[order(average_foreign_pop_world_regions$mean, decreasing = T),]

require(scales)
ggplot(subset(immi_by_year_by_world_region, 
              variable %in% average_foreign_pop_world_regions$variable[1:7]), 
       aes(x=year, y=foreign_residents, colour=variable)) +
  geom_line() +
  geom_point() + 
  geom_line(data=subset(immi_by_year_by_world_region, 
                        variable %in% average_foreign_pop_world_regions$variable[6:22]), 
            aes(x=year, y=foreign_residents, group=variable), alpha=0.4) + 
  geom_label(data=subset(immi_by_year_by_world_region, 
                         variable %in% average_foreign_pop_world_regions$variable[1:7] &
                           year == 2016), 
             aes(x=year, y=foreign_residents, colour=variable, label = variable),
             alpha = 0.5, nudge_x = 1.2) +
  scale_x_continuous(limits = c(2003, 2018)) + 
  scale_y_continuous(label = scales::comma) +
  theme_bw() +
  labs(x=NULL) +
  guides(colour=FALSE)


```

```{r, cache=T}
## FEMALE

istat_immi_prov_fem <-
  read.csv("istat_provincial_foreign_residents.csv", 
           comment.char="#", encoding = "latin1", stringsAsFactors = FALSE)
istat_immi_prov_fem <- 
  subset(istat_immi_prov_fem, Sesso == "femmine")
istat_immi_prov_fem[["Tipo.di.indicatore.demografico"]] <- NULL
istat_immi_prov_fem[["Sesso"]] <- NULL
istat_immi_prov_fem[["Flags"]] <- NULL
names(istat_immi_prov_fem)[4] <- 'Residenti'

# Drop macro-regions
macro_regions <- c('Nord-est', 'Nord-ovest', 'Centro', 'Sud', 'Italia',
                   'Isole')
istat_immi_prov_fem <-
  subset(istat_immi_prov_fem, !(Territorio %in% macro_regions))

# Drop regions
codici_regioni <- read.csv("~/ownCloud/gabriele_lega_immi/data/istat_codice_regione/codici_regioni.csv")
regioni <- as.character(codici_regioni$denominazione_regione)
regioni <- gsub("/(.*)", "", regioni)
regioni <- regioni[-which(regioni == "VALLE D'AOSTA")]

regioni[!tolower(regioni) %in% tolower(istat_immi_prov_fem$Territorio)]

istat_immi_prov_fem <- 
  subset(istat_immi_prov_fem, !(tolower(Territorio) %in% tolower(regioni)))

# Drop double
istat_immi_prov <- 
  subset(istat_immi_prov_fem, !(grepl('Provincia Autonoma', 
                                      Territorio)))

istat_immi_prov_fem$Territorio <- 
  gsub(" /(.*)", "", istat_immi_prov_fem$Territorio)

require(countrycode)
edited_nomi_provenienza <-
  read.csv("name_foreign_countries.csv",
           stringsAsFactors = FALSE)
dict_names <- 
  edited_nomi_provenienza$En
names(dict_names) <- edited_nomi_provenienza$It

istat_immi_prov_fem$country_machine_ready <- 
  istat_immi_prov_fem$Paese.di.cittadinanza

istat_immi_prov_fem$country_machine_ready <- 
  dict_names[istat_immi_prov_fem$country_machine_ready]

istat_immi_prov_fem$iso3c <-
  countrycode(istat_immi_prov_fem$Paese.di.cittadinanza,
              "country.name", "iso3c")

istat_immi_prov_fem[is.na(istat_immi_prov_fem$iso3c),]$iso3c <-
  countrycode(istat_immi_prov_fem[is.na(istat_immi_prov_fem$iso3c),]$country_machine_ready,
              "country.name", "iso3c")

# unique(istat_immi_prov_fem[is.na(istat_immi_prov_fem$iso3c),]$Paese.di.cittadinanza)

istat_immi_prov_fem[istat_immi_prov_fem$Paese.di.cittadinanza=="Kosovo",]$iso3c <- "KSV"
istat_immi_prov_fem[istat_immi_prov_fem$Paese.di.cittadinanza=="Serbia e Montenegro",]$iso3c <-"SCG"

istat_immi_prov_fem$world_region <- 
  countrycode(istat_immi_prov_fem$iso3c, 'iso3c', 'region')

istat_immi_prov_fem[istat_immi_prov_fem$Paese.di.cittadinanza=="Kosovo",]$world_region <- "Southern Europe"
istat_immi_prov_fem[istat_immi_prov_fem$Paese.di.cittadinanza=="Serbia e Montenegro",]$world_region <-"Southern Europe"
istat_immi_prov_fem[istat_immi_prov_fem$Paese.di.cittadinanza=="Taiwan (ex Formosa)",]$world_region <- "Eastern Asia"

istat_immi_prov_fem <-
  istat_immi_prov_fem %>%
  dplyr::group_by(Anno, world_region) %>%
  dplyr::summarize(foreign_residents = sum(Residenti))

istat_immi_prov <-
  read.csv('istat_provincial_immi_prov_wt_musl_and_world_region_pop.csv')

# Male/female ratio
require(dplyr)
tot_immigration <-
  istat_immi_prov %>%
  dplyr::group_by(year) %>%
  dplyr::summarize(Eastern.Europe = sum(Eastern.Europe),
                   Southern.Europe = sum(Southern.Europe),
                   Northern.Africa = sum(Northern.Africa),
                   South.America = sum(South.America),
                   Southern.Asia = sum(Southern.Asia))

require(reshape2)
tot_immigration <- melt(tot_immigration, id.vars = 'year')
tot_immigration$variable <- as.character(tot_immigration$variable)

names(istat_immi_prov_fem) <- c('year', 'variable', 'value')
istat_immi_prov_fem$variable <- gsub(" ", "\\.", istat_immi_prov_fem$variable)

tot_immigration <- merge(tot_immigration, istat_immi_prov_fem, by = c('year','variable'),
                         all = FALSE)

tot_immigration$ratio <- 
  with(tot_immigration, value.y / value.x)

ggplot(tot_immigration, aes(x=year, y=ratio, colour=variable)) +
  geom_line() +
  geom_point() + 
  geom_label(data=subset(tot_immigration, year == 2016), 
             aes(x=year, y=ratio, colour=variable, label = variable),
             alpha = 0.5, nudge_x = 1.2) +
  scale_x_continuous(limits = c(2003, 2018)) + 
  scale_y_continuous(label = scales::percent) +
  theme_bw() +
  labs(x=NULL, y="Percentage of female immigrants") +
  guides(colour=FALSE)

```



```{r}
# Maps
require("rgdal")
require("maptools")
require(plyr)
require(dplyr)

shape <- readOGR(".",
                 layer="istat_Prov2008_WGS84_g_nord_simp")
shape@data$id <- rownames(shape@data)
gpclibPermit()
shape.points <- fortify(shape, region="id")
shape.df <- join(shape.points, shape@data, by="id")

# Calculate distance from centroids
# trueCentroids <- gCentroid(shape, byid=TRUE)
# point_df <- as.data.frame(trueCentroids@coords)
# point_df$group <- 1
#   
# ggplot(shape.df) + 
#   aes(long, lat, group=group) + 
#   geom_polygon(fill=NA) +
#   geom_path(color="black", size = 0.1) +
#   coord_equal() +
#   geom_point(data = point_df, aes(x=x,y=y, group=group))
# require(geosphere)
# distance_of_centroids <- 
#   dist2Line(
#     spTransform(trueCentroids[1], CRS("+proj=longlat +datum=WGS84")),
#     spTransform(shape, CRS("+proj=longlat +datum=WGS84")))
# mean(distance_of_centroids[,1])
# [1] 12791.94

shape_entire <- 
  readOGR(".",
          layer="istat_Prov2008_WGS84_g")

require(rgeos)
neigh <- gTouches(shape_entire, byid=TRUE)
colnames(neigh) <- shape_entire$COD_PRO
rownames(neigh) <- shape_entire$COD_PRO
diag(neigh) <- TRUE


# Add neigh average 
# immi_votes$neigh_foreign_residents <- NA
# immi_votes$neigh_muslim_foreign_residents <- NA
# immi_votes$neigh_Numero.votanti <- NA
# immi_votes$neigh_Voti.LN <- NA

for (i in which(immi_votes$code_region < 9)) {
  y <- immi_votes$year[i]
  p <- immi_votes$code_province[i]
  neigh_provinces <- colnames(neigh)[neigh[as.character(p),]]
  
  tmp_df <- subset(immi_votes, year == y & code_province %in% neigh_provinces)
  
  immi_votes$neigh_foreign_residents[i] <- sum(tmp_df$foreign_residents)
  immi_votes$neigh_female_foreign_residents[i] <- sum(tmp_df$female_foreign_residents)
  immi_votes$neigh_male_foreign_residents[i] <- sum(tmp_df$male_foreign_residents)
  
  immi_votes$neigh_muslim_foreign_residents[i] <- sum(tmp_df$muslim_foreign_residents)
  immi_votes$neigh_Numero.votanti[i] <- sum(tmp_df$Numero.votanti)
  immi_votes$neigh_Voti.LN[i] <- sum(tmp_df$Voti.LN)
  
  for (r in region_vars) {
    immi_votes[[paste0('neigh_', r)]][i] <-
      sum(tmp_df[[r]])
  }
  
}

immi_votes$neigh_perc_lega <- 
  with(immi_votes, neigh_Voti.LN / neigh_Numero.votanti)
immi_votes$neigh_perc_pop_straniera <- 
  with(immi_votes, neigh_foreign_residents / neigh_Numero.votanti)
immi_votes$neigh_perc_pop_female <- 
  with(immi_votes, neigh_female_foreign_residents / neigh_Numero.votanti)
immi_votes$neigh_perc_pop_male <- 
  with(immi_votes, neigh_male_foreign_residents / neigh_Numero.votanti)
immi_votes$neigh_perc_pop_musul <- 
  with(immi_votes, neigh_muslim_foreign_residents / neigh_Numero.votanti)
immi_votes$neigh_perc_pop_non_musul <- 
  with(immi_votes, neigh_foreign_residents - neigh_muslim_foreign_residents / neigh_Numero.votanti)

for (r in region_vars) {
  immi_votes[[paste0('neigh_perc_', r)]] <-
    immi_votes[[paste0('neigh_', r)]] / immi_votes$neigh_Numero.votanti
}

```


```{r}
yrs <-c(2004, 2005, 2006, 2008, 2009, 2010, 2013, 2014, 2015)

require("RColorBrewer")
require("viridis")
require(scales)
lega_pal <- colorRampPalette(brewer.pal(9, "YlGn"))
sc_lega <- 
  scale_fill_gradientn(colours = lega_pal(100), 
                       limits=c(0, max(immi_votes$perc_lega, na.rm=T)),
                       na.value = "white", labels=percent)
immi_pal <- colorRampPalette(viridis(15))
sc_immi <- 
  scale_fill_gradientn(colours = rev(viridis(100)), 
                       limits=c(0, max(immi_votes$perc_pop_straniera, na.rm=T)),
                       na.value = 'white', labels=percent)
musl_pal <- colorRampPalette(magma(15))
sc_musl <- 
  scale_fill_gradientn(colours = rev(magma(100)), 
                       limits=c(0, max(immi_votes$perc_pop_musul, na.rm=T)),
                       na.value = 'white', labels=percent)

musl_over_foreign_pal <- colorRampPalette(magma(15))
sc_musl_over_foreign <- 
  scale_fill_gradientn(colours = rev(magma(100)), 
                       limits=c(0, max(immi_votes$perc_pop_musul / 
                                         immi_votes$perc_pop_straniera, na.rm=T)),
                       na.value = 'white', labels=percent)

unempl_pal <- colorRampPalette(plasma(15))
sc_unempl <- 
  scale_fill_gradientn(colours = rev(plasma(100)), 
                       limits=c(0, 
                                max(subset(immi_votes,
                                           code_region < 9)$perc_unempl/100, na.rm=T)),
                       na.value = 'white', labels=percent)

extractLegend <- function(a.gplot) {
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

map_theme <-
  theme(axis.line=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        # panel.background=element_rect(fill = 'white'),
        panel.border=element_blank(),
        panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        plot.background=element_blank(),
        legend.margin=unit(0, "cm"),
        legend.position="bottom")

map_lega <- list()
map_immi <- list()
map_musl <- list()
map_musl_over_foreign <- list()
map_unempl <- list()

for (y in yrs) {
  tmp_shp_df <- 
    merge(shape.df, subset(immi_votes, year == y), by.x = 'COD_PRO', by.y = 'code_province')
  
  map_lega[[as.character(y)]] <-
    ggplot(tmp_shp_df) + 
    aes(long, lat, group=group, fill=perc_lega) + 
    geom_polygon() +
    geom_path(color="black", size = 0.1) +
    coord_equal() +
    theme_bw() +
    sc_lega +
    labs(title = y, fill = "Lega Nord votes (% of voters)") +
    map_theme +
    theme(plot.margin=unit(c(-1,0,-1,0),"cm"))
  
  legend_lega <- extractLegend(map_lega[[as.character(y)]])
  
  map_lega[[as.character(y)]] <-
    ggplot_gtable(ggplot_build(map_lega[[as.character(y)]] + theme(legend.position="none")))
  
  map_immi[[as.character(y)]] <-
    ggplot(tmp_shp_df) + 
    aes(long, lat, group=group, fill=perc_pop_straniera) + 
    geom_polygon() +
    geom_path(color="black", size = 0.1) +
    coord_equal() +
    theme_bw() +
    sc_immi +
    labs(title = y, fill = "Extra-EU residents (% of electorate)") +
    map_theme +
    theme(plot.margin=unit(c(-1,0,-1,0),"cm"))
  
  legend_immi <- extractLegend(map_immi[[as.character(y)]])
  
  map_immi[[as.character(y)]] <-
    ggplot_gtable(ggplot_build(map_immi[[as.character(y)]] + theme(legend.position="none")))
  
  map_musl[[as.character(y)]] <-
    ggplot(tmp_shp_df) + 
    aes(long, lat, group=group, fill=perc_pop_musul) + 
    geom_polygon() +
    geom_path(color="black", size = 0.1) +
    coord_equal() +
    theme_bw() +
    sc_musl +
    labs(title = y, fill = "Extra-EU Muslims (% of electorate)") +
    map_theme +
    theme(plot.margin=unit(c(-1,0,-1,0),"cm"))
  
  legend_musl <- extractLegend(map_musl[[as.character(y)]])
  
  map_musl[[as.character(y)]] <-
    ggplot_gtable(ggplot_build(map_musl[[as.character(y)]] + theme(legend.position="none")))
  
  map_musl_over_foreign[[as.character(y)]] <-
    ggplot(tmp_shp_df) + 
    aes(long, lat, group=group, fill=perc_pop_musul / perc_pop_straniera) + 
    geom_polygon() +
    geom_path(color="black", size = 0.1) +
    coord_equal() +
    theme_bw() +
    sc_musl_over_foreign +
    labs(title = y, fill = "Extra-EU Muslims (% of foreign population)") +
    map_theme +
    theme(plot.margin=unit(c(-1,0,-1,0),"cm"))
  
  legend_musl_over_foreign <- extractLegend(map_musl_over_foreign[[as.character(y)]])
  
  map_musl_over_foreign[[as.character(y)]] <-
    ggplot_gtable(ggplot_build(map_musl_over_foreign[[as.character(y)]] + theme(legend.position="none")))
  
  map_unempl[[as.character(y)]] <-
    ggplot(tmp_shp_df) + 
    aes(long, lat, group=group, fill=perc_unempl/100) + 
    geom_polygon() +
    geom_path(color="black", size = 0.1) +
    coord_equal() +
    theme_bw() +
    sc_unempl +
    labs(title = y, fill = "Unemployment (% of population over 15)") +
    map_theme +
    theme(plot.margin=unit(c(-1,0,-1,0),"cm"))
  
  legend_unempl <- extractLegend(map_unempl[[as.character(y)]])
  
  map_unempl[[as.character(y)]] <-
    ggplot_gtable(ggplot_build(map_unempl[[as.character(y)]] + theme(legend.position="none")))
  
}

require(gridExtra)
grob1 <- arrangeGrob(map_lega[[1]],map_lega[[2]],map_lega[[3]],map_lega[[4]],map_lega[[5]],
                     map_lega[[6]], map_lega[[7]], map_lega[[8]], map_lega[[9]], nrow=1)
grob2 <- arrangeGrob(map_immi[[1]],map_immi[[2]],map_immi[[3]],map_immi[[4]],map_immi[[5]],
                     map_immi[[6]], map_immi[[7]], map_immi[[8]], map_immi[[9]], nrow=1)
grob3 <- arrangeGrob(map_musl[[1]],map_musl[[2]],map_musl[[3]], map_musl[[4]],map_musl[[5]],
                     map_musl[[6]], map_musl[[7]], map_musl[[8]], map_musl[[9]], nrow=1)
grob4 <- arrangeGrob(map_musl_over_foreign[[1]],map_musl_over_foreign[[2]],map_musl_over_foreign[[3]], map_musl_over_foreign[[4]],map_musl_over_foreign[[5]],
                     map_musl_over_foreign[[6]], map_musl_over_foreign[[7]], map_musl_over_foreign[[8]], map_musl_over_foreign[[9]], nrow=1)
grob5 <- arrangeGrob(map_unempl[[1]],map_unempl[[2]],map_unempl[[3]], map_unempl[[4]],map_unempl[[5]],
                     map_unempl[[6]], map_unempl[[7]], map_unempl[[8]], map_unempl[[9]], nrow=1)

row1 <- arrangeGrob(grob1, nrow=1)
row2 <- arrangeGrob(grob2, nrow=1)
row3 <- arrangeGrob(grob3, nrow=1)
row4 <- arrangeGrob(grob4, nrow=1)
row5 <- arrangeGrob(grob5, nrow=1)
```

```{r, fig.width = 20, fig.height = 12.5}
       grid.arrange(row1, legend_lega, row2, legend_immi, row3, legend_musl, row4, legend_musl_over_foreign, row5, legend_unempl, 
                    nrow = 10,  heights = unit(c(3, 2, 4, 2, 4, 2, 4, 2, 4, 1), "cm"))
```


```{r, results = 'asis'}
# Regression

require(plyr)

# Rename variables
immi_votes <- plyr::rename(immi_votes, 
                     c("perc_lega" = "votes",
                       "Numero.elettori" = "voters",
                       "affluenza"="tornout",
                       "perc_unempl"="unemp",
                       "perc_lega_reg"="regionalVotes",
                       "neigh_perc_lega"="neighVotes",
                       "perc_pop_straniera"="foreignPop",
                       "perc_pop_female"="femaleForeignPop",
                       "perc_pop_male"="maleForeignPop",
                       "neigh_perc_pop_straniera"="neighForeignPop",
                       "neigh_perc_pop_female"="neighFemaleForeignPop",
                       "neigh_perc_pop_male"="neighMaleForeignPop",
                       "perc_pop_musul"="muslimForeignPop",
                       "neigh_perc_pop_musul"="neighMuslimForeignPop",
                       "perc_pop_non_musul"="nonMuslimForeignPop",
                       "neigh_perc_pop_non_musul"="neighNonMuslimForeignPop",
                       "perc_Eastern.Europe"="easternEuropePop",
                       "perc_Southern.Europe"="southernEuropePop",
                       "perc_Northern.Africa"="northernAfricaPop",
                       "perc_South.America"="southAmericaPop",
                       "perc_Southern.Asia"="southernAsiaPop",
                       "neigh_perc_Eastern.Europe"="neighEasternEuropePop",
                       "neigh_perc_Southern.Europe"="neighSouthernEuropePop",
                       "neigh_perc_Northern.Africa"="neighNorthernAfricaPop",
                       "neigh_perc_South.America"="neighSouthAmericaPop",
                       "neigh_perc_Southern.Asia"="neighSouthernAsiaPop"))

immi_votes$voters <-  immi_votes$voters/1000
immi_votes$unemp <- immi_votes$unemp/100

require(plm)
fe.mod1 <- 
  plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp,
      index = c("year","code_province"), model = "within", data = immi_votes)
fe.mod2 <- 
  plm(votes ~ 
        voters + tornout + 
        neighVotes +
        unemp,
      index = c("year","code_province"), model = "within", data = immi_votes)
fe.mod3 <- 
  plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        neighForeignPop,
      index = c("year","code_province"), model = "within", data = immi_votes)
fe.mod4 <- 
  plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        neighNonMuslimForeignPop +
        neighMuslimForeignPop,
      index = c("year","code_province"), model = "within", data = immi_votes)
fe.mod5 <- 
  plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        over65Pop + 
        neighNonMuslimForeignPop +
        neighMuslimForeignPop,
      index = c("year","code_province"), model = "within", data = immi_votes)

pool.mod5 <- 
  plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        over65Pop + 
        neighNonMuslimForeignPop +
        neighMuslimForeignPop,
      index = c("year","code_province"), model = "pooling", data = immi_votes)

# Test for multicollinearity
require(car)
vif(lm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        over65Pop + 
        neighNonMuslimForeignPop +
        neighMuslimForeignPop,
       data = immi_votes))

fe.mod6 <- 
    plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        over65Pop +
        neighMaleForeignPop +
        neighFemaleForeignPop,
      index = c("year","code_province"), model = "within", data = immi_votes)

fe.mod7 <- 
    plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        neighFemaleForeignPop +
        neighMuslimForeignPop,
      index = c("year","code_province"), model = "within", data = immi_votes)

fe.mod8 <- 
    plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        neighFemaleForeignPop +
        neighMuslimForeignPop +
        over65Pop,
      index = c("year","code_province"), model = "within", data = immi_votes)

fe.mod9 <-
  plm(votes ~ 
        voters + tornout + 
        regionalVotes +
        unemp +
        over65Pop + 
        neighForeignPop +
        neighEasternEuropePop +
        neighSouthernEuropePop +
        neighNorthernAfricaPop +
        neighSouthAmericaPop +
        neighSouthernAsiaPop,
      index = c("year","code_province"), model = "within", data = immi_votes)



require(stargazer)
stargazer(fe.mod1, fe.mod2, fe.mod3, fe.mod4, fe.mod5,
          font.size = 'footnotesize', omit.stat = 'f')
stargazer(fe.mod6, fe.mod7, fe.mod8, fe.mod9,
          font.size = 'footnotesize', omit.stat = 'f')
stargazer(fe.mod5, pool.mod5,
          font.size = 'footnotesize', omit.stat = 'f')
```


```{r}
require(car)
scatterplot(femaleForeignPop / foreignPop ~ over65Pop, data = immi_votes)
```

```{r}
scatterplotMatrix(
   ~ votes +
    tornout + regionalVotes + unemp +
    over65Pop + neighNonMuslimForeignPop + neighMuslimForeignPop,
  data = immi_votes, cex = 0.5)
```

